#pragma once
#include "Image.hpp"
#include "ImageHelper.hpp"
#include <Function/GaussianFunction.hpp>
#include <Function/Gaussian2DFunction.hpp>
#include <Math/DVector.hpp>
#include <Math/DMatrix.hpp>

namespace zzz{
template<typename T>
class ImageFilter
{
public:
//general convolve
  static void ConvolveRow(DVectord & kern, Image<T> * src, Image<T> * dst);  //convolve row
  static void ConvolveCol(DVectord & kern, Image<T> * src, Image<T> * dst);  //convolve column
  static void ConvolveDouble(DVectord &kern, Image<T> *src, Image<T> * dst);  //first convolve row, then convolve column
  static void Convolve2D(DMatrixd &kern, Image<T> *src, Image<T> *dst);    //convolve a 2d block
private:
  static T ConvolveRowHelper(DVectord & kernel, Image<T> * src, zuint r, zuint c);
  static T ConvolveColHelper(DVectord & kernel, Image<T> * src, zuint r, zuint c);
  static T Convolve2DHelper(DMatrixd & kernel, Image<T> * src, zuint r, zuint c);

public:
  //Gaussian Blur
  static void GaussBlurImage(Image<T> * src, Image<T> * dst, int kernelsize, double sigma);
  static DVectord GaussianKernel1D(int kernelsize, double sigma);
  static DMatrixd GaussianKernel2D(int kernelsize, double sigma);


};
//general
template<typename T>
T ImageFilter<T>::ConvolveRowHelper(DVectord & kernel, Image<T> * src, zuint r, zuint c) 
{
  T pixel(0);
  int cen = kernel.size / 2;
  for (int i = 0; i < (int)kernel.size; i++) 
  {
    int col = c + (i - cen);
    if (col < 0) col = 0;
    if (col >= (int)src->Cols()) col = src->Cols() - 1;
    pixel += src->At(r,col)*kernel[i];
  }
  return PixelHelper<T>::Clamp(pixel);
}

template<typename T>
void ImageFilter<T>::ConvolveRow(DVectord & kern, Image<T> * src, Image<T> * dst) 
{
  dst->SetSize(src->Size());
  for (zuint r = 0; r < src->Rows(); r++) for (zuint c = 0; c < src->Cols(); c++) 
    dst->At(r, c)=ConvolveRowHelper(kern, src, r, c);
}

template<typename T>
T ImageFilter<T>::ConvolveColHelper(DVectord & kernel, Image<T> * src, zuint r, zuint c)
{
  T pixel(0);
  int cen = kernel.size / 2;
  for (int j = 0; j < (int)kernel.size; j++) 
  {
    int row = r + (j - cen);
    if (row < 0) row = 0;
    if (row >= (int)src->Rows()) row = src->Rows() - 1;
    pixel += src->At(row,c)*kernel[j];
  }
  return PixelHelper<T>::Clamp(pixel);
}

template<typename T>
void ImageFilter<T>::ConvolveCol(DVectord & kern, Image<T> * src, Image<T> * dst) 
{
  dst->SetSize(src->Size());
  for (zuint r = 0; r < src->Rows(); r++) for (zuint c = 0; c < src->Cols(); c++) 
    dst->At(r, c)=ConvolveColHelper(kern, src, r, c);
}


template<typename T>
void ImageFilter<T>::ConvolveDouble(DVectord & kern, Image<T> * src, Image<T> * dst) 
{
  dst->SetSize(src->Size());
  Image<T> tmpImage(src->Size());
  ConvolveRow(kern, src, &tmpImage);
  ConvolveCol(kern, &tmpImage, dst);
}


template<typename T>
T zzz::ImageFilter<T>::Convolve2DHelper(DMatrixd & kernel, Image<T> * src, zuint r, zuint c)
{
  T pixel(0);
  int cenr = kernel.nrow / 2;
  int cenc = kernel.ncol / 2;
  for (int i=0; i<(int)kernel.nrow; i++) 
  {
    int row=r+(i-cenr);
    if (row<0) row=0;
    if (row>=(int)src->Rows()) row=src->Rows()-1;
    for (int j=0; j<(int)kernel.ncol; j++)
    {
      int col=c+(j-cenc);
      if (col<0) col=0;
      if (col>=(int)src->Cols()) col=src->Cols()-1;
      pixel+=src->At(row,col)*kernel(i,j);
    }
  }
  return PixelHelper<T>::Clamp(pixel);  
}

template<typename T>
void zzz::ImageFilter<T>::Convolve2D(DMatrixd &kern, Image<T> *src, Image<T> *dst)
{
  dst->SetSize(src->Size());
  for (zuint r=0; r<src->Rows(); r++) for (zuint c=0; c<src->Cols(); c++)
    dst->At(r,c)=Convolve2DHelper(kern,src,r,c);
}
//////////////////////////////////////////////////////////////////////////
//for gaussian
template<typename T>
void ImageFilter<T>::GaussBlurImage(Image<T> * src, Image<T> * dst, int kernelsize, double sigma) 
{
//These two method should product the same result
//method1: do 2 times 1d convolve(row and column)
  DVectord convkernel = GaussianKernel1D(kernelsize, sigma);
  ImageFilter<T>::ConvolveDouble(convkernel,src,dst);
//method2: do a 2d convolve
//  DMatrixd convkernel=GaussianKernel2D(kernelsize, sigma);
//  ImageFilter<T>::Convolve2D(convkernel,src,dst);
}
template<typename T>
DVectord ImageFilter<T>::GaussianKernel1D(int kernelsize, double sigma) 
{
  GaussianFunction func(0,sigma);
  if (kernelsize%2==0) kernelsize++;
  DVectord kern(kernelsize);
  int c = kernelsize / 2;
  for (int i = 0; i < (kernelsize + 1) / 2; i++) 
  {
    double v = func(i);
    kern[c+i] = v;
    kern[c-i] = v;
  }
  kern/=kern.Sum();
  return kern;
}

template<typename T>
DMatrixd ImageFilter<T>::GaussianKernel2D(int kernelsize, double sigma) 
{
  Gaussian2DFunction func(Vector2d(0,0),sigma);
  if (kernelsize % 2 == 0) kernelsize++;
  DMatrixd mat(kernelsize,kernelsize);
  int c = kernelsize / 2;
  for (int i = 0; i < (kernelsize + 1) / 2; i++) for (int j = 0; j < (kernelsize + 1) / 2; j++) 
  {
    double v = func(Vector2d(i,j));
    mat(c+i,c+j) = v;
    mat(c-i,c+j) = v;
    mat(c+i,c-j) = v;
    mat(c-i,c-j) = v;
  }
  mat/=mat.Sum();
  return mat;
}

//////////////////////////////////////////////////////////////////////////


}
